Competing Nuclear Quantum Effects and Hydrogen-Bond Jumps in Hydrated Kaolinite

Recent work has shown that the dynamics of hydrogen bonds in pure clays are affected by nuclear quantum fluctuations, with different effects for the hydrogen bonds holding different layers of the clay together and for those within the same layer. At the clay–water interface there is an even wider range of types of hydrogen bond, suggesting that the quantum effects may be yet more varied. We apply classical and thermostated ring polymer molecular dynamics simulations to show that nuclear quantum effects accelerate hydrogen-bond dynamics to varying degrees. By interpreting the results in terms of the extended jump model of hydrogen-bond switching, we can understand the origins of these effects in terms of changes in the quantum kinetic energy of hydrogen atoms during an exchange. We also show that the extended jump mechanism is applicable not only to the hydrogen bonds involving water, but also those internal to the clay.

* sı Supporting Information ABSTRACT: Recent work has shown that the dynamics of hydrogen bonds in pure clays are affected by nuclear quantum fluctuations, with different effects for the hydrogen bonds holding different layers of the clay together and for those within the same layer. At the clay−water interface there is an even wider range of types of hydrogen bond, suggesting that the quantum effects may be yet more varied. We apply classical and thermostated ring polymer molecular dynamics simulations to show that nuclear quantum effects accelerate hydrogen-bond dynamics to varying degrees. By interpreting the results in terms of the extended jump model of hydrogen-bond switching, we can understand the origins of these effects in terms of changes in the quantum kinetic energy of hydrogen atoms during an exchange. We also show that the extended jump mechanism is applicable not only to the hydrogen bonds involving water, but also those internal to the clay.
H ydrogen bonding is central to the properties of clays, 1−4 with H-bonds not only holding together the layers of clay materials but also controlling their interactions with surrounding solvents such as water. Clay−water systems show a rich variety of hydrogen-bond types and strengths. These interactions are crucial in processes such as the adsorption of solute on the surface of the clay, 5,6 its wetting, 7−9 and the swelling of the material by intercalation between the layers. 10,11 Recent work by the authors 4 has shown that the properties of dry clays can be affected by nuclear quantum effects (NQEs) such as zero-point energy, leading to weaker hydrogen bonds. The magnitude of these quantum effects depends on whether the H-bonds held together different clay layers or acted within a single layer; in clay−water systems, there is a wider range of H-bonding, with clays accepting H-bonds from the surrounding water, and in some cases able to donate hydrogen bonds to the water. The range of H-bond strengths exhibited by clay− water systems 1−3 raises the possibility that NQEs may affect different types of hydrogen bond to different degrees.
In this Letter, we use classical and ring polymer molecular dynamics to show that NQEs do indeed affect the hydrogen bonds in clay−water systems to different degrees. By interpreting our results using the angular jump model of hydrogen-bond switching, 12−14 we are able to understand these quantum effects in terms of the change in confinement of a water molecule as it undergoes a jump. We show further that the angular jump mechanism is present in H-bonds that do not involve water, hinting at its universality.
We focus on NQEs in hydrated kaolinite [Al 2 Si 2 O 5 (OH) 4 ], a 1:1-type dioctahedral clay whose layers comprise an octahedral sheet containing aluminum and a tetrahedral sheet containing silicon. Pure kaolinite contains H-bonds both within a clay layer and between neighboring layers; the kaolinite−water interface adds H-bonds donated by water molecules to the silica sheets and H-bonds between water and the aluminol sheets, in which either water or clay O−H groups may be the donor. 2,3 As in ref 4, since our main focus is on the importance of NQEs, we use the CLAYFF-TRPMD force field developed in that work, which extends the popular CLAYFF force field 15−18 to give accurate vibrational spectra in path integral simulations. Similarly, the q-TIP4P/F force field was used to model water; 19 this model accounts well for the experimental properties of pure water in path integral simulations, without requiring a significant computational expense. Classical molecular dynamics (MD), path-integral molecular dynamics (PIMD), 20 and thermostated ring polymer molecular dynamics (TRPMD) 21 simulations were carried out using the i-PI 22,23 and LAMMPS 24 codes, with NQEs accounted for by the PIMD and TRPMD calculations. To avoid the problem of artificial electric fields at kaolinite−water interfaces described in ref 25, two types of system were simulated: one in which both clay layers were oriented with the aluminol sheets toward the water and one with both silica sheets oriented toward water. Example simulation inputs as well as further details are given in the Supporting Information.
To identify events in which an H-bond donor changed its acceptor, the stable-states picture of chemical reactions was used, 26 as in refs 12 and 27. That is, a hydrogen-bond switch occurs whenever a donor oxygen O, whose initial acceptor is O a , forms a hydrogen bond with a new acceptor O b . We use a strict geometrical criterion for H-bonding, as described in the Supporting Information. The different types of H-bond studied in this letter are given in Figure 1. We began by calculating the hydrogen-bond lifetime correlation functions, 27 with the average over all possible H-bonds. n R (t) is 1 if the Hbond is intact at time t and 0 otherwise, and n P (t) is 1 if the original H-bond is no longer intact at time t, but another Hbond has been formed with the same donor O−H group. Absorbing boundary conditions are used, so that once a hydrogen bond is broken and a new one formed with the same donor, n R (t) = 0 thereafter. By fitting to an exponential decay, C(t) ≃ e −t/τ 0 , we find the hydrogen-bond exchange time τ 0 . In path integral simulations the n R (t) and n P (t) functions are evaluated using the ring-polymer centroid. Figure 2 shows the lifetime correlation function C(t) from MD and TRPMD calculations for the various types of H-bond that are present in these simulations. In Figure 2b, the dynamics of intralayer hydrogen bonds are split into those of H-bonds that remain in the minimum-energy structure (i.e., at zero temperature, denoted T = 0 K) and those of H-bonds where the donor would participate in interlayer H-bonds at zero temperature but at finite temperatures is participating in an intralayer H-bond (denoted T > 0 K). This distinction was shown to be important in ref 4, in which NQEs on the two types of intralayer hydrogen bonds were shown to be very different. The resulting H-bond exchange times τ 0 are shown in Table 1 and are in accordance with existing results which show that hydrogen bonding between pairs of water molecules is weaker than with the (hydrophilic) aluminol surface and stronger than with the (hydrophobic) silica surface. 1,2 The exchange times increase in the order water−silica < water− water < water−aluminol. Table 1 shows that hydrogen bonds donated by aluminol O−H groups and accepted by water are weaker than those donated by water and accepted by the aluminol layer, with the latter being 50% longer lived. Due to the higher electronegativity of Al than of H, the partial negative charge of oxygen in water is greater than that in an aluminol O−H group, making water the stronger H-bond donor. This difference in hydrogen-bonding strength was previously observed in ab initio MD simulations. 2,3 For intralayer Hbonds, there are two distinct time scales: T > 0 K H-bonds are much shorter lived than intralayer H-bonds that persist at zero temperature. For the T = 0 K intralayer H-bonds, the activation of an exchange involves the rearrangement of the relatively rigid solid structure, meaning that its characteristic time is the longest observed. In Figure 2, some correlation functions are not single exponentials, particularly for water−  The Journal of Physical Chemistry Letters pubs.acs.org/JPCL Letter aluminol and finite-temperature intralayer H-bonds. This behavior has been observed previously for hydrogen-bonding dynamics at water−mineral interfaces, and attributed to a distribution of exchange times τ 0 . 28,29 Although all of the quantum effects for H-bond exchange times are small, there is a variation among the different types. For the H-bonds involving water, we note that the longer the classical exchange time τ CL , the larger the quantum effect, with water−silica H-bonds having the smallest value for each of these and water−aluminol H-bonds the largest. This is in accord with previous results suggesting that more strongly bound systems, with longer lifetimes, are more susceptible to NQEs, in the context of ion solvation 30,31 and proton transfer reactions. 32 Since clays are a significant source of contaminant removal from soil, 6 this result suggests that NQEs may increase the rate of this process. However, the results for kaolinite's internal H-bonds do not fit into this trend; although the T > 0 K intralayer H-bonds are longer-lived than water−silica Hbonds, there are essentially no NQEs within the error bars. In addition, despite having the same classical exchange time, water−water and interlayer H-bonds incur different quantum effects. Crucially, the differing ratios of classical to quantum exchange times implies that NQEs in clay−water systems cannot be straightforwardly accounted for by running classical simulations at elevated temperatures.
To better understand the different magnitudes of quantum effect seen for hydrated kaolinite we note that the time τ 0 is a key ingredient in the extended jump model of H-bond switching, in which hydrogen bonds change their acceptors through a large-angle jump of the donor O−H group; 14 this has been observed in a range of systems, 12,33−35 including mineral−water interfaces. 36−38 While it has been shown that water molecules donating hydrogen bonds to a material also change their acceptor via a large-angle jump, it is not yet known whether hydrogen bonds internal to clay minerals also undergo these exchanges, allowing us to investigate the applicability of this mechanism to systems beyond water. For H-bonds involving water, the type is formatted "Donor− acceptor". Results are shown for classical molecular dynamics (labelled "CL") and thermostated ring polymer molecular dynamics (labelled "NQE"). The ratio of the two is also shown. b τ CL and τ NQE are the same within their error bars. While the trajectories of H-bonds involving water look qualitatively very similar to those observed in pure water, 12,27 the donor−acceptor distances inside the clay (panels (d)−(f) in Figure 3) behave nonmonotonically, and the behavior of This asymmetry reflects the fact that 95% of interlayer H-bonds become intralayer H-bonds after undergoing a jump, and 25% of intralayer H-bonds become interlayer H-bonds. Jumps where both the initial and final acceptor are in the same layer are much more symmetric, as shown in the Supporting Information. This discrepancy may appear at first to be at odds with detailed balance, indicating that jumps do not happen independently�rather, the equilibrium distribution of Hbonds is maintained by reversals of jumps. The populations of interlayer and intralayer H-bonds in dry kaolinite, according to our previous work, 4 are 0.51:0.19:0.10:0.20 interlayer:T = 0 K intralayer:T > 0 K intralayer:dangling (i.e., without an H-bond acceptor). The nonmonotonic behavior is less straightforward to understand: future work will focus on investigating whether this is due to the complexity of the mechanism or to limited statistics. The ratios of jumps to different types of H-bond acceptor are also instructive: for water molecules initially donating an H-bond to the aluminol surface, the ratio of final acceptors that are on the aluminol surface to those that are water molecules is around 1:1; for donation to the silica surface, this ratio is closer to 1:3, illustrating the higher strength of H-bonds formed with the aluminol. As shown in the Supporting Information, NQEs make very little difference to the trajectories of H-bond jumps including water, in accord with the results of ref 39. To understand the origin of these NQEs, we computed the centroid virial kinetic energy tensor,   Table 1. The inset illustrates this tensor as an ellipsoid during a jump event: at the transition state, the H atom is able to move toward either the initial or final acceptor, meaning its confinement is low in the perpendicular direction; it is only weakly H-bonded, meaning that it is relatively confined in the parallel direction.
The Journal of Physical Chemistry Letters pubs.acs.org/JPCL Letter acceptors is relatively low, while decreases, indicating a greater delocalization in this direction as the O−H bond can move just as easily toward the initial and the final acceptor.
The inset of Figure 4 shows this in more detail. Qualitatively, these results are in accord with previous work. 4 These two effects combine to give the total quantum kinetic energy, which generally decreases at the transition state. Table 2 shows the changes in kinetic energy on going to the transition state, decomposed into parallel and perpendicular components. does not account for the quantum effect observed by itself, as the average potential energy will also change, but a greater decrease in the reaction barrier due to kinetic energy correlates fairly well with larger quantum effects in τ 0 (see Table 1).
A significant cause of the quantum effects in H-bond exchange times is the change in confinement of the H atom, leading to a decrease in the free energy barrier to undergoing a jump. The effects of motion parallel and perpendicular to the O−H bond compete with each other, with the latter dominating in all cases. This competition of effects has been observed before for hydrogen-bonded systems, 19,31,40 with the dominant effect being determined by the type of motion that contributes most to H-bond breaking. 40 This means that, while the NQEs in our simulations are extremely subtle, they act as a direct probe of the change in geometry that occurs on going to the transition state. Since only a few hundred jump events were collected for intralayer H-bonds that are only present at finite temperatures, the kinetic energy trajectories are extremely noisy (see the Supporting Information), meaning that the apparent lack of quantum effects in their exchange time τ 0 cannot yet be fully understood. As in Figure 3, the symmetry in the quantum kinetic energy trajectories indicates the types of H-bonds that are preferentially formed during a jump: Hbonds initially donated from water to the silica layer have highly asymmetric trajectories because half of the final acceptors are water molecules; for this reason, the section of the water−Si trajectory with t > 0 is very similar to that of the water−water trajectory. On the other hand, H-bonds donated to the aluminol layer have much more symmetric trajectories, since a greater proportion of the final acceptor is also in this layer.
In this Letter, we have shown that hydrated clays exhibit a variety of nuclear quantum effects, which are rationalized using the angular jump mechanism of hydrogen-bond exchange. The different degrees of quantum effects are found to be a manifestation of the change in confinement that an H-bond experiences at the jump transition state. Future work will focus on understanding the complex jump mechanism in and around clays further, paying attention to the nonmonotonicity in trajectories, as well as the effect of quantum fluctuations at the interface between materials and ionic solutions, 41 on surfaces at which water molecules dissociate, 42 and on using more sophisticated potentials 43